﻿//===================================================================
//
//            本代码由David Wang原创，仅用于学习之用
//        All these codes can only be used to learning 
//
//          Perlin噪声1D、2D生成算法，向 Perlin 致敬
//                       2017.08.15                
//
//===================================================================
using System;
using System.Collections.Generic;

namespace DavidWang
{
    class PerlinNoise1D
    {
        static double pointA;
        static double pointB;
        public  float Amplitude { set; get; }
        public  float Frequency { set; get; }
        private static int Wavelength = 1;
        PRNG myPRNG = new PRNG();
        public PerlinNoise1D()
        {
            pointA = myPRNG.Next();
            pointB = myPRNG.Next();
            Amplitude = 200.0f;
            Frequency = 0.005f;
            Wavelength = 1;
        }
        private double Cos_Interpolate(double a, double b, double t)
        {
            double ft = t * Math.PI;
            t = (1 - Math.Cos(ft)) * 0.5;
            return a * (1 - t) + t * b;
        }
        public int Perlin(int x)
        {
            Wavelength = (int)(1.0f / Frequency);
            double y;
            if (x % Wavelength == 0)
            {
                pointA = pointB;
                pointB = myPRNG.Next();
                y = pointA * Amplitude;
            }
            else
            {
                y = Cos_Interpolate(pointA, pointB, (double)(x % Wavelength) / Wavelength) * Amplitude;
            }
            return (int)y;
        }
    }

    class PRNG
    {
        static long M = 4294967296;
        static long A = 1664525;
        static long C = 1;
        static long Z;
        public PRNG()
        {
            Random rdm = new Random();
            Z = (long)Math.Floor((rdm.NextDouble() * M));
        }
        public double Next()
        {
            Z = (A * Z + C) % M;
            return (double)Z / (double)M;
        }
    }

    class PerlinNoise2D
    {
        #region 常用变量
        public float Amplitude { set; get; }
        public float Frequency { set; get; }
        private static float[,] gradients;
        private static readonly int[] permutation = { 151,160,137,91,90,15,               
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,    
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
    };
        private static  int[] p;
        //perlin 噪声函数中使用的变量都定义在这里，以提高效率
        private int Xi, Yi, grad11, grad12, grad21, grad22;
        private float fracX, fracY,
              noise11, noise12, noise21, noise22,
              interpolatedx1, interpolatedx2;
        #endregion
        public PerlinNoise2D()
        {
            Amplitude = 1.0f;
            Frequency = 1.0f;
            gradients = new float[8,2];
            GenerateGradients();
            p = new int[256];
            Random rdm = new Random();
            int j;
            for (int x = 0; x < 256; x++)
            {
                j = (int)rdm.Next(256);
                p[x] = permutation[j];
            }
            
            
        }

        public float Perlin(int x, int y)
        {

            Xi = GetFloor((float)x * Frequency);
            Yi = GetFloor((float)y * Frequency);
            fracX = (float)x * Frequency - (float)Xi;
            fracY = (float)y * Frequency - (float)Yi;
            grad11 = p[(Xi + p[Yi & 255]) & 255] & 7;
            grad12 = p[(Xi + 1 + p[Yi & 255]) & 255] & 7;
            grad21 = p[(Xi + p[(Yi + 1) & 255]) & 255] & 7;
            grad22 = p[(Xi + 1 + p[(Yi + 1) & 255]) & 255] & 7;
            noise11 = DotProduct(gradients[grad11,0], gradients[grad11, 1], fracX, fracY);
            noise12 = DotProduct(gradients[grad12,0], gradients[grad12, 1], fracX - 1.0f, fracY);
            noise21 = DotProduct(gradients[grad21,0], gradients[grad21, 1], fracX, fracY - 1.0f);
            noise22 = DotProduct(gradients[grad22,0], gradients[grad22, 1], fracX - 1.0f, fracY - 1.0f);
            fracX = Fade(fracX);
            fracY = Fade(fracY);
            interpolatedx1 = (float)Interpolate(noise11, noise12, fracX);
            interpolatedx2 = (float)Interpolate(noise21, noise22, fracX);
            
            return (float)Interpolate(interpolatedx1, interpolatedx2, fracY) * Amplitude;
        }
        #region Helper
        private void GenerateGradients()
        {
            for (int i = 0; i < 8; ++i)
            {
                gradients[i,0] = (float)Math.Cos(0.785398163f * (float)i);  // ( 2 * PI / 8) * i
                gradients[i,1] = (float)Math.Sin(0.785398163f * (float)i); 
            }
        }        
        private int GetFloor(float value)
        {
            return (value >= 0 ? (int)value : (int)value - 1);
        }
        float DotProduct(float gradX,float gradY,float distanceX, float distanceY)
        {
            return (gradX * distanceX + gradY * distanceY);
        }
        private double Interpolate(double a, double b, double t)
        {
            return a * (1 - t) + t * b;
        }
        private float Fade(float x)
        {            
            return (x * x * x * (x * (6 * x - 15) + 10));  //缓和函数 6x^5 - 15x^4 + 10x^3 
        }
        #endregion
    }

    class PerlinNoise3D
    {
        #region 常用变量
        public float Amplitude { set; get; }
        public float Frequency { set; get; }
        private static List<Vector3> gradients;
        private static readonly int[] permutation = { 151,160,137,91,90,15,
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
    };
        private static int[] p;
        //perlin 噪声函数中使用的变量都定义在这里，以提高效率
        private int Xi, Yi,Zi, aaa, aba, aab, abb, baa, bba, bab, bbb;
        private float fracX, fracY, fracZ, u, v, w;
        double x1, x2, y1, y2;
        #endregion
        public PerlinNoise3D()
        {
            Amplitude = 1.0f;
            Frequency = 1.0f;
            GenerateGradients();
            p = new int[512];
            for (int x = 0; x < 512; x++)
            {
                p[x] = permutation[x % 256];
            }
        }
        public float Perlin(int x, int y, int z)
        {
            Xi = GetFloor((float)x * Frequency);
            Yi = GetFloor((float)y * Frequency);
            Zi = GetFloor((float)z * Frequency);
            fracX = (float)x * Frequency - (float)Xi;
            fracY = (float)y * Frequency - (float)Yi;
            fracZ = (float)z * Frequency - (float)Zi;
            u = Fade(fracX);
            v = Fade(fracY);
            w = Fade(fracZ);
            aaa = p[p[p[Xi]     + Yi]      + Zi];
            aba = p[p[p[Xi]     + Yi + 1]  + Zi];
            aab = p[p[p[Xi]     + Yi]      + Zi + 1];
            abb = p[p[p[Xi]     + Yi + 1]  + Zi + 1];
            baa = p[p[p[Xi + 1] + Yi]      + Zi];
            bba = p[p[p[Xi + 1] + Yi + 1]  + Zi];
            bab = p[p[p[Xi + 1] + Yi]      + Zi + 1];
            bbb = p[p[p[Xi + 1] + Yi + 1]  + Zi + 1];
 
            x1 = Interpolate(DotProduct(gradients[aaa & 0xF], fracX, fracY, fracZ), DotProduct(gradients[baa & 0xF], fracX - 1, fracY, fracZ), u);
            x2 = Interpolate(DotProduct(gradients[aba & 0xF], fracX, fracY - 1, fracZ), DotProduct(gradients[bba & 0xF], fracX - 1, fracY - 1, fracZ), u);
            y1 = Interpolate(x1, x2, v);
            x1 = Interpolate(DotProduct(gradients[aab & 0xF], fracX, fracY, fracZ - 1), DotProduct(gradients[bab & 0xF], fracX - 1, fracY, fracZ - 1), u);
            x2 = Interpolate(DotProduct(gradients[abb & 0xF], fracX, fracY - 1, fracZ - 1), DotProduct(gradients[bbb & 0xF], fracX - 1, fracY - 1, fracZ - 1), u);
            y2 = Interpolate(x1, x2, v);

            return (float)Interpolate(y1, y2, w) * Amplitude;
        }

        #region Helper
        private void GenerateGradients()
        {
            gradients = new List<Vector3>
            {
                new Vector3(1,1,0),
                new Vector3(-1,1,0),
                new Vector3(1,-1,0),
                new Vector3(-1,-1,0),
                new Vector3(1,0,1),
                new Vector3(-1,0,1),
                new Vector3(1,0,-1),
                new Vector3(-1,0,-1),
                new Vector3(0,1,1),
                new Vector3(0,-1,1),
                new Vector3(0,1,-1),
                new Vector3(0,-1,-1),            
                
                new Vector3(1,1,0),
                new Vector3(0,-1,1),
                new Vector3(-1,1,0),
                new Vector3(0,-1,-1)
            };
        }
        class Vector3
        {
            public int x { get; }
            public int y { get; }
            public int z { get; }
            public Vector3(int xi, int yi, int zi)
            {
                x = xi;
                y = yi;
                z = zi;
            }
        }

        private double DotProduct(Vector3 v3, float fracX, float fracY, float fracZ)
        {
            return v3.x * fracX + v3.y * fracY + v3.z * fracZ;
        }
        private int GetFloor(float value)
        {
            return (value >= 0 ? (int)value : (int)value - 1);
        }
        private double Interpolate(double a, double b, double t)
        {
            return a * (1 - t) + t * b;
        }
        private float Fade(float x)
        {
            return (x * x * x * (x * (6 * x - 15) + 10));  //缓和函数 6x^5 - 15x^4 + 10x^3 
        }
        #endregion
    }

    class PerlinNoise4D
    {
        #region 常用变量
        public float Amplitude { set; get; }
        public float Frequency { set; get; }
        private static List<Vector4> gradients;
        private static readonly int[] permutation = { 151,160,137,91,90,15,
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
    };
        private static int[] p;
        //perlin 噪声函数中使用的变量都定义在这里，以提高效率
        private int Xi, Yi, Zi,Wi, aaaa, aaba, aaab, aabb, abaa, abba, abab, abbb, baaa, baba, baab, babb, bbaa, bbba, bbab, bbbb;
        private float fracX, fracY, fracZ, fracW, r, g, b, a;
        double x1, x2, y1, y2, y3, y4, z1, z2;
        #endregion
        
        public PerlinNoise4D()
        {
            Amplitude = 1.0f;
            Frequency = 1.0f;
            GenerateGradients();
            p = new int[512];
            for (int x = 0; x < 512; x++)
            {
                p[x] = permutation[x % 256];
            }
        }
        public float Perlin(int x, int y, int z,int w)
        {
            Xi = GetFloor((float)x * Frequency);
            Yi = GetFloor((float)y * Frequency);
            Zi = GetFloor((float)z * Frequency);
            Wi = GetFloor((float)w * Frequency);
            fracX = (float)x * Frequency - (float)Xi;
            fracY = (float)y * Frequency - (float)Yi;
            fracZ = (float)z * Frequency - (float)Zi;
            fracW = (float)w * Frequency - (float)Wi;
            r = Fade(fracX);
            g = Fade(fracY);
            b = Fade(fracZ);
            a = Fade(fracW);

            aaaa = p[p[p[p[Xi] + Yi]     + Zi]     + Wi];
            aaba = p[p[p[p[Xi] + Yi]     + Zi + 1] + Wi];
            aaab = p[p[p[p[Xi] + Yi]     + Zi]     + Wi + 1];
            aabb = p[p[p[p[Xi] + Yi]     + Zi + 1] + Wi + 1];
            abaa = p[p[p[p[Xi] + Yi + 1] + Zi]     + Wi];
            abba = p[p[p[p[Xi] + Yi + 1] + Zi + 1] + Wi];
            abab = p[p[p[p[Xi] + Yi + 1] + Zi]     + Wi + 1];
            abbb = p[p[p[p[Xi] + Yi + 1] + Zi + 1] + Wi + 1];

            baaa = p[p[p[p[Xi + 1] + Yi]     + Zi]     + Wi];
            baba = p[p[p[p[Xi + 1] + Yi]     + Zi + 1] + Wi];
            baab = p[p[p[p[Xi + 1] + Yi]     + Zi]     + Wi + 1];
            babb = p[p[p[p[Xi + 1] + Yi]     + Zi + 1] + Wi + 1];
            bbaa = p[p[p[p[Xi + 1] + Yi + 1] + Zi]     + Wi];
            bbba = p[p[p[p[Xi + 1] + Yi + 1] + Zi + 1] + Wi];
            bbab = p[p[p[p[Xi + 1] + Yi + 1] + Zi]     + Wi + 1];
            bbbb = p[p[p[p[Xi + 1] + Yi + 1] + Zi + 1] + Wi + 1];

            x1 = Interpolate(DotProduct(gradients[aaaa & 0x1F], fracX, fracY, fracZ, fracW), DotProduct(gradients[baaa & 0x1F], fracX - 1, fracY, fracZ, fracW), r);
            x2 = Interpolate(DotProduct(gradients[aaba & 0x1F], fracX, fracY, fracZ-1,fracW), DotProduct(gradients[baba & 0x1F], fracX - 1, fracY, fracZ -1 ,fracW), r);
            y1 = Interpolate(x1, x2, g);
            x1 = Interpolate(DotProduct(gradients[aaab & 0x1F], fracX, fracY, fracZ,fracW -1), DotProduct(gradients[baab & 0x1F], fracX - 1, fracY, fracZ,fracW -1), r);
            x2 = Interpolate(DotProduct(gradients[aabb & 0x1F], fracX, fracY, fracZ - 1,fracW -1), DotProduct(gradients[babb & 0x1F], fracX - 1, fracY, fracZ - 1, fracW -1), r);
            y2 = Interpolate(x1, x2, g);

            x1 = Interpolate(DotProduct(gradients[abaa & 0x1F], fracX, fracY-1, fracZ, fracW), DotProduct(gradients[bbaa & 0x1F], fracX - 1, fracY-1, fracZ, fracW), r);
            x2 = Interpolate(DotProduct(gradients[abba & 0x1F], fracX, fracY-1, fracZ - 1, fracW), DotProduct(gradients[bbba & 0x1F], fracX - 1, fracY-1, fracZ - 1, fracW), r);
            y3 = Interpolate(x1, x2, g);
            x1 = Interpolate(DotProduct(gradients[abab & 0x1F], fracX, fracY-1, fracZ, fracW - 1), DotProduct(gradients[bbab & 0x1F], fracX - 1, fracY-1, fracZ, fracW - 1), r);
            x2 = Interpolate(DotProduct(gradients[abbb & 0x1F], fracX, fracY-1, fracZ - 1, fracW - 1), DotProduct(gradients[bbbb & 0x1F], fracX - 1, fracY-1, fracZ - 1, fracW - 1), r);
            y4 = Interpolate(x1, x2, g);

            z1 = Interpolate(y1, y2, b);
            z2 = Interpolate(y3, y4, b);

            return (float)Interpolate(z1, z2, a) * Amplitude;
        }
        #region Helper
        private void GenerateGradients()
        {
            gradients = new List<Vector4>
            {
                new Vector4(0,-1,-1,-1),
                new Vector4(0,-1,-1,1),
                new Vector4(0,-1,1,-1),
                new Vector4(0,-1,1,1),
                new Vector4(0,1,-1,-1),
                new Vector4(0,1,-1,1),
                new Vector4(0,1,1,-1),
                new Vector4(0,1,1,1),

                new Vector4(-1,-1,0,-1),
                new Vector4(-1,1,0,-1),
                new Vector4(1,-1,0,-1),
                new Vector4(1,1,0,-1),
                new Vector4(-1,-1,0,1),
                new Vector4(-1,1,0,1),
                new Vector4(1,-1,0,1),
                new Vector4(1,1,0,1),

                new Vector4(-1,0,-1,-1),
                new Vector4(1,0,-1,-1),
                new Vector4(-1,0,-1,1),
                new Vector4(1,0,-1,1),
                new Vector4(-1,0,1,-1),
                new Vector4(1,0,1,-1),
                new Vector4(-1,0,1,1),
                new Vector4(1,0,1,1),

                new Vector4(0,-1,-1,0),
                new Vector4(0,-1,-1,0),
                new Vector4(0,-1,1,0),
                new Vector4(0,-1,1,0),
                new Vector4(0,1,-1,0),
                new Vector4(0,1,-1,0),
                new Vector4(0,1,1,0),
                new Vector4(0,1,1,0)

            };
        }
        class Vector4
        {
            public int x { get; }
            public int y { get; }
            public int z { get; }
            public int w { get; }
            public Vector4(int xi, int yi, int zi,int wi)
            {
                x = xi;
                y = yi;
                z = zi;
                w = wi;
            }
        }

        private double DotProduct(Vector4 v4, float fracX, float fracY, float fracZ, float fracW)
        {
            return v4.x * fracX + v4.y * fracY + v4.z * fracZ + v4.w * fracW;
        }
        private int GetFloor(float value)
        {
            return (value >= 0 ? (int)value : (int)value - 1);
        }
        private double Interpolate(double a, double b, double t)
        {
            return a * (1 - t) + t * b;
        }
        private float Fade(float x)
        {
            return (x * x * x * (x * (6 * x - 15) + 10));  //缓和函数 6x^5 - 15x^4 + 10x^3 
        }
        #endregion
    }
}
